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Abstract. We propose a variational framework for the resolution of a non¬ 
hydrostatic Saint-Venant type model with bottom topography. This model is 
a shallow water type approximation of the incompressible Euler system with 
free surface and slightly differs from the Green-Nagdhi model, see [12] for more 
details about the model derivation. 

The numerical approximation relies on a prediction-correction type scheme 
initially introduced by Chorin-Temam [16] to treat the incompressibility in the 
Navier-Stokes equations. The hyperbolic part of the system is approximated 
using a kinetic finite volume solver and the correction step implies to solve a 
mixed problem where the velocity and the pressure are defined in compatible 
finite element spaces. 

The resolution of the incompressibility constraint leads to an elliptic prob¬ 
lem involving the non-hydrostatic part of the pressure. This step uses a varia¬ 
tional formulation of a shallow water version of the incompressibility condition. 

Several numerical experiments are performed to confirm the relevance of 
our approach. 


1. Introduction. Starting from the incompressible Euler or Navier-Stokes system, 
the hydrostatic assumption consists in neglecting the vertical acceleration of the 
fluid. More precisely - and with obvious notations - the momentum along the 
vertical axis of the Euler equation 


dw dw 

■ 

reduces in the hydrostatic context to 


dw dp 


= -g, 


dp 

dz 


= -g- 


( 1 ) 


Such an assumption produces important consequences over the structure and com¬ 
plexity of the model. Indeed, Eq. (1) implies that the pressure p is no longer the 
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Lagrange multiplier of the incompressibility constraint and p can be expressed, for 
free surface flows, as a function of the water depth of the fluid. Therefore, the 
hydrostatic assumption implies that the resulting model, even though it describes 
an incompressible fluid, has common features with models arising in compressible 
fluid mechanics. 

In geophysical problems, the hydrostatic assumption coupled with a shallow wa¬ 
ter type description of the flow is often used. Unfortunately, these models do not 
represent phenomena containing dispersive effects for which the non-hydrostatic 
contribution cannot be neglected. And more complex models have to be consid¬ 
ered to take into account this kind of phenomena, together with numerical methods 
able to discretize the high order derivative terms coming from the dispersive ef¬ 
fects. Many shallow water type dispersive models have been proposed such as 
KdV, Boussinesq, Green-Naghdi, see [21, 14, 6, 30, 31, 27, 19, 26, 2, 3, 13]. The 
modeling of the non-hydrostatic effects for shallow water flows does not raise insu¬ 
perable difficulties but their discretization is more tricky. Numerical techniques for 
the approximaion of these models have been recently proposed [15, 11, 28]. 

The model studied in the present paper has been derived and studied in [12]. Its 
numerical approximation based on a projection-correction strategy [16] is described 
in [1] . In [1] , the discretization of the elliptic part arising from the non-hydrostatic 
terms is carried out in a finite difference framework. It is worth noticing that the 
numerical scheme given in [1] is endowed with robustness and stability properties 
such as positivity, well-balancing, discrete entropy and wet/dry interfaces treatment. 

The main contents of this paper is the derivation and validation of the correction 
step in a variational framework. Since the derivation in a 2d context of the model 
proposed in [12] does not raise difficulty, the results depicted in this paper pave the 
way for a discretization of the 2d model over an unstructured mesh, and we will 
often maintain general notations as far as possible. 

Notice that the non-hydrostatic model we consider slightly differs from the well- 
known Green-Naghdi model [21] but the numerical approximation proposed in this 
paper can also be used for the numerical approximation of the Green-Naghdi system. 

Let U C M, be a Id domain (an interval) and T = T^n U Tout its boundary (see 
figure 1). The non-hydrostatic model derived in [12, 1] reads 


dHu 


dH dHu 

dt dx 

= 0, 

( 2 ) 

d 

+ 3^ + HPnh) 

= -{gH + 2p„f^)^, 

( 3 ) 

dHw dHwu 

= ‘2Pnh^ 

( 4 ) 

dt dx 

dHu d{H + 2zb) 

U -\-2w 

ox ox 

= 0, 

( 5 ) 


where H is the water depth, zt the topography and Pnh the non-hydrostatic part of 
the pressure. The variables denoted with a bar recall that this model is obtained 
performing an average along the water depth of the incompressible Euler system 
with free surface. The velocity field is denoted u = {u,wY with u (resp. w) the 
horizontal (resp. vertical) component. 
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We denote r] = H + Zf, the free surface of the fluid. In addition, we give the 
following notations 


n = 


n 

0 


( 6 ) 


with n the unit outward normal vector at T (in Id, n = ± 1), n represents the unit 
outward normal vector of the domain covered by the fluid, namely x [z},,r]]. We 
also consider the gradient operator 


Vo = 



( 7 ) 


Figure 1. Notations and domain definition. 


The smooth solutions of the system (2)-(5) satisfy moreover an energy conserva¬ 
tion law 

with 

-_ H{u^+ w^) ^ gH{r] + Zb) ^ 

Note that (5) represents a shallow water version of the divergence free constraint, 
for which the non hydrostatic pressure pnt plays the role of a Lagrange multiplier. 
Notice that considering pnh = 0 and neglecting (4), the system (2)-(3),(5) reduces 
to the classical Saint-Venant system. 

The paper is organized as follows. First we give a rewriting of the model and 
we present the prediction-correction method, the main part being the variationnal 
formulation of the correction part. Then in Section 3, we detail the numerical 
approximation. Finally, in Section 4, numerical simulations validating the proposed 
discretization techniques are presented. 


2. The projection scheme for the non-hydrostatic model. Projection meth¬ 
ods have been introduced by A. Chorin and R. Temam [33] in order to compute 
the pressure for incompressible Navier-Stokes equations. These methods, based on 
a time splitting scheme, have been widely studied and applied to treat the incom¬ 
pressibility constraint (see [24, 35, 34]). We develop below an analogue of this 
method for shallow water flow. In order to describe the fractional time step method 
we use, we propose a rewritting of the model (2)-(5). 


2 . 1 . 


A rewritting. Let us introduce the two operators Vg. 

■d{H+2zb)\ 


and divg, 


^SW f - 


divg^, (v) = 


dx ^ j dx 


dHvx 


- 2 / 
d{H 


vi- 


)■ 

2zb) 


+ 2112 , 


defined by 

( 10 ) 

( 11 ) 


dx dx 

with V = {vi,V 2 Y- We assume for a while that / and v are smooth enough. The 
shallow water form of the divergence operator div^^j (resp. of the gradient operator 
Vgu,) corresponds to a depth averaged version of the divergence (resp. gradient) 
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appearing in the incompressible Euler and Navier-Stokes equations. Notice that 
the two operators Vsm , divs^, defined by (lO)-(ll) are H and Zb dependent and we 
assume that H and zt are sufficiently smooth functions. One can check that these 
operators verify the fundamental duality relation 


/ divsu, (v)/ dx = - \7sw f ■'vdx+ [Hvi /]p . 
J ^ J Q, 


These definitions allow to rewrite the model (2)-(5) as 

dH dHu 


dx 


= 0 , 


dHu d 


dt + ^ = -gHWoZb, 

divs^ (u) = 0, 


( 12 ) 


(13) 

(14) 

(15) 


with Vo defined by (7). 

The system (13)-(15) can be written in the compact form 


where we denote 


and 


X = 


Rnh — 


dX 

Ih 


= S{X), 

(16) 


divsu, (u) 

= 0, 

(17) 

H 

\ / 

Hu \ 


Hu 

, nx) = 

Hu^ -H |i72 1 

(18) 

Hw 

/ V 

Huw ) 



0 

^sw Pnh 


S{X) = 


0 


(19) 


-gHXoZb 

Let be given time steps At" and note t" = detailed in [1], the 

projection scheme for system (16)-(17) consists in the following time splitting 


j^n+1/2 


Xn+l 


divsu, u"+^ 


A" - At" ^FiX'^) + At" ^(A"), 

OX 

A"+i/2_ At" i?0+\ 


0 , 


( 20 ) 

( 21 ) 

( 22 ) 


with u"+i = ) . 

The first two equations of (20) consist in the classical Saint-Venant system with 
topography and the third equation is an advection equation for the quantity Hw. 
Equations (21)-(22) describe the correction step allowing to determine the non hy¬ 
drostatic part of the pressure and hence giving the corrected state A"+^. The 
numerical resolution of (20) - especially the first two equations - has received an 
extensive coverage and efficient and robust numerical techniques exist, mainly based 
on finite volume approach, see [8, 5]. The derivation of a robust and efficient nu¬ 
merical technique for the resolution of the correction step (21)-(22) is the key point. 
A strategy based on a finite difference approach has been proposed, studied and 
validated in [1]. Unfortunately, the finite difference framework does not allow to 
tackle situations with unstructured meshes in 2 or 3 dimensions. It is the key point 
of this paper to propose a variational formulation of the correction step coupled 
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with a finite volume discretization of the prediction step. 

The time discretization in the numerical scheme described above corresponds 
to a fractional time step strategy with a first order Euler scheme, explicit for the 
hyperbolic part and implicit for the elliptic part. For hyperbolic conservation laws, 
the second-order accuracy in time is usually recovered by the Heun method [7, 
8] that is a slight modification of the second order Runge-Kutta method. More 
precisely, for a dynamical system written under the form 


1 = 

(23) 

the Heun scheme consists in dehning by 


n 1 ~n+2 

n +1 _y +y 
y 2 

(24) 

with 


y"+i = y" + Arw), 

yn+2 ^ ~n+l ArF(l)”+l). 


The model we have to discretize has the general form 


f 

(25) 

B{X) = 0, 


and we propose the numerical scheme 


2 

(26) 

with the two steps defined by 


= X{r) + ArF(A",:p”+i), 

(27) 

H(X”+i) = 0, 

(28) 

and 


j^n+2 ^ j^n+1 ArE(A"+\p”+2), 

(29) 

H(A"+2) = 0, 

(30) 


where and are the solutions of the elliptic equations derived from the 
divergence free constraints (28) and (30) respectively. Notice that (27) is a compact 
form of the fractional scheme (20)-(21) where the intermediate step no more appears. 

It is easy to prove the scheme (26) is second order accurate. Indeed assuming F 
and B are enough smooth, we have 

X{e + AF) = X(F) + AFX{F) + ^ ^ ^ X{e) F 0{{AFf), 

or equivalently using (25) 

X{F + At") = X{F) + ArF{X{F),p{F)) 

2 \ oX op at J ^ ' 

(31) 
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Now from (27) and (29), we get 


^n+2 


xn ArF(2s:”,:p”+i) + ArF(X"+\:p"+2) 
X" + 2ArF(A",:p"+i) 


+ (Ar)2 

+o{iArf). 




9F(A”,p"+i) 


9F(A”,p"+i) dp^+^ 
dp dt 


( 32 ) 


Using (26), we see that relations (31) and (32) are equivalent up to third order 
terms. 


2.2. The correction step. In this part, we consider we have at our disposal a 
space discretization of Eq. (20) solving the hydrostatic part of the model and we 
focus on the correction step (21)-(22). 


2.2.1. Variational formulation. The correction step (21)-(22) writes, 


H-n+i 

(i7u)"+i+ArV«„p((+i 

div^u, (u”+^) 


_^n+l/2, 

(i7u)"+i/^ 

0 . 


(33) 

(34) 

(35) 


For the sake of clarity, in the following we will drop the notation with a bar and 
we denote p instead of Pnh- Likewise we drop the superscript for the corrected 
states. 

Equations (34)-(35) is a mixed problem in velocity/pressure, its approximation 
leads to the variational mixed problem 
Find p G Q, u e Vo with 

Q = {ge L2(fl)|V«^gG (L2(fl))2}, (36) 


and 


Vb = {v = {vi,V 2 ) e (L^(f2))^|divs^ (v) G L^{Vl),vi\y = 0}, 
such that 

f (iLu + At"Vsu,p) • vdcc = f (iLu)”'*'^^^ • V dec, VvgVq, 

Jq Jq 

/ divs„(u)gdx = 0, '^qeQ. 

Jn 

Introducing the bilinear forms 

a(u,v) = / Hu vdx, Vu, vGVb, 

Jn 

b{v,q) = - diVswMqdx, VvGUo,VqGQ, 

Jn 

the problem (38)-(39) becomes: 

Find p G Q and u G Vb such that 

^a(u,v)+6(v,p) = ^a(u"+i/^v), Vv G Vq, 

b{u,q) = 0, yqGQ. 


(37) 

(38) 

(39) 


(40) 

(41) 
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2.2.2. The pressure equation. Formally, we take v in the form v = with q G Qo 
and 


Qo 

Then, with (12) and (41), we get 

a(u,v) 


= {q€ Qqjr = 0 }. 


/ u ■ Vsw qdx 

Jq 


= - diVsu, (u) • qdx 

Jn 


= 0 . 


If we introduce the shallow water version of the Laplacian operator Asw 


AswP = divsu, 



it is natural to consider the new variational formulation: 
Find VL GV,p G Qo,sw such that 


^a(u,v) + (v, Vs^p) = ^a(u"+^/^v) Vv e 1/, (42) 

{Asn,P,q) = ^ (divs„ (u”+i/2),g) , WqGQo.sw, (43) 

with Qo,sw = {<7 G QoldiVsu, G T^(f^)}. 

From (43), we deduce 

^sw (p) = ^dWsw (u"+^/2^. (44) 

The resolution of the equation (44) allows to update the velocity at the correction 
step (34). 

Notice that the equation (44) is equivalent to apply the operator divg^j to the 
equation ( 21 ) and to use the shallow water free divergence ( 22 ) to eliminate u. 

Remark 1. Taking the functional spaces V and Qo^sw, the problem (42)-(43) is 
not equivalent to the problem (40)-(41). 

Remark 2. Notice that Eq. (44) has the form of a Sturm-Liouville type equation. 


2.2.3. Inf-sup condition. To ensure that the saddle problem (40)-(41) is well posed, 
the Babuska-Brezzi [9, 32] condition 


37 > 0 , 7 < inf 


&(v,g) 


has to be satisfied. Denoting by B the weak operator defined by 'dv G Vq, Bv = 
b{v,q),yq G Q, we have 


kerS* = {qGQ\B^q^0} 


= {9 G Q| / y/sw g • V dx = 0 Vv e Do}- 

Jn 

We assume H > > 0. Because of the positivity of id, it is obvious that the 

bilinear form a is coercive. Choosing v G (L^(D))^ and q such as Vgu, q G (L^(r2))^, 
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it follows that Wsw 9 = 0, then q = 0 and ker = 0. Indeed, in contrast with Navier- 
Stokes equations for which the pressure is defined up to an additive constant, the 
non hydrostatic of the shallow water equations is fully defined. Therefore, the mixed 
problem (40)-(41) satisfies the inf-sup condition and admits a unique solution. 

2.3. Boundary conditions. In this section, we still consider that the hydrostatic 
part is provided and we study the compatibility of the boundary conditions be¬ 
tween the hydrostatic part and the projection part. Therefore, the compatibility 
between the pressure and velocity at boundary needs to be studied. To this aim, 
we first provide the conditions required to impose Dirichlet or Neumann pressure 
at boundary, and then, we couple these conditions with the hydrostatic part. 

We consider a more general case taking the space V 


V = {v = {vi,V2) e (L^(fl))^|div^,^ V e L^(n)}, 
and we introduce the bilinear form 

c(v,p) = J Hpvinds, \/vgV,p€Q, 

with n the unit outward normal vector defined by (6). In one dimension, c{'v,p) = 


Therefore instead of (40)-(41), we consider the problem: 
Find u e y, p G Q such that. 


1 

At" 


a(u,v) -I- &(v,p) 
6(u, q) 


^a(u”+i/^v)-hc(v,p), VvgF, 
0, \/q G Q. 


(45) 

(46) 


Notice that divs^, (u) = Vq ■ [H^) + u- (Ug -f Ub) and VswP = H'^o{p) — p(ns + nb) 
with Vo defined by (7) and Ug (resp. Ub) the (non-unit) normal vector at the 
surface (resp. at the bottom) 



Moreover, we have 


nb 


dzh 

dx 

1 


/ divsu, (u) dx 

Jq 


Hunds 


u • (ug -f Ub) dx. 


Hence, to satisfy the divergence free condition, the velocity u should verify 


Hunds 


/ u • (Ug -f Ub) dx. 

Jn 


Dirichlet condition for the pressure. From the variational formulation (40)-(41) of 
the projection scheme, a natural boundary condition for the pressure is a Dirichlet 
condition. At ri,(z = in, out), pjc, = Po then c{v,p) = fj. povin ds and we take 
V £ V, q e Qi with 


Qt = 


{q G L\n)\Vsyz q G (L"(D))^ (;|r. = 0}. 
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Neumann boundary for the pressure. The Neumann boundary condition for the 
projection scheme is not natural and to enforce such a condition, the elliptic prob¬ 
lem (42)-(43) has to be considered. Taking now q G Qsnn with Q^w = {q G 
Qldivsu, (-^ 75 ^) G the problem is rewritten 

^a(u,v)-b (v,Vs^p) = ^a(u"+i/^v) Vv G 

{A,^p,q) = ^ (|diw^ -b ^5(u,g), '^qeQsw, 

with c, the bilinear form 


c(u,p) = J {Hu + At"‘\7 swP\i) qndj — J qnd'y. 

Many studies have been done to choose an appropriate variational formulation for 
this problem. In [23] J-L. Guermond explores the different variational formulations 
in order to enforce a Neumann pressure boundary condition, in [25] some equiv¬ 
alent formulations are given to switch between Neumann and Dirichlet boundary 
conditions. 

Taking the normal component at the boundary T^ of the momentum equation, 
it follows that 


I I I .OJtL . . ^ ( |^+l/2 I \ 


,dH 


H 


dn 


dn 


At" 


We note §f-|ri = A, * = in, out. 

• In case /3i = 0, a Neumann boundary condition for the pressure is deduced of 
a Dirichlet condition for u. 

dp I 


op. i , |n-H/2 I S 

^\ t . = ip:(“lr, -»lr.). 

• In the other cases, it gives a mixed boundary condition 


(47) 


fjlr.+A p = (48) 

Then, in the two cases, we have imposed a Dirichlet velocity condition, that leads 
to take V & Vi and q G Q, with for i = in, out 

Vi = {-v = {vi,V 2 ) G (L^(fl))^|divsu, (v) G L^(f2),i;i|r, = 0}. (49) 

Let us now give the coupling boundary conditions between the prediction step and 
the correction step. Indeed, in the projection part, boundary conditions need to be 
set in order to be consistent with the hydrostatic part. 

Concerning the prediction step, we consider the well known Saint-Venant system 
and we assume that the Riemann invariant remains constant along the associated 
characteristic. This approach has been introduced in [10] and distinguishes fluvial 
and torrential boundaries depending on the Froude number Fr = —. Usual bound¬ 
ary conditions consist to impose a flux qo at the inflow boundary and a water depth 
at the outflow boundary. It is also classical to let a free outflow boundary, setting 
a Neumann boundary condition for the water depth and for the velocity. In both 
cases, we give the boundary conditions that have to be set in the correction step. 
We consider the first situation in which we set a flux at the inflow T^^ and a given 
depth at the outflow Tout- Assuming a fluvial flow, this case consists in solving 
a Riemann problem at the interface Tin where the global flux is given by qo = 
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(< 701 , 902 )* = That gives the boundary values Hq = 

uo = HTi /2 ^0 = !1'+i/2 from the hyperbolic part. This leads to obtain a 

Dirichlet condition for the pressure at the left boundary of the correction part. 
Moreover, if H is given for the outflow, we preconize to give a mixed condition for 
the pressure that corresponds to the boundary condition (48) 


P\i 


dp 

dn 


dH 

+ 

dn 


= 0 , 

= 0 , 


that leads to take u G Vout, with the definition (49) and 

P&Qin = {qG L^{n)\V,^qG L^{n), 9|r.„ = 0}. 

We now consider the second situation in which we still impose a flux qo at the 
inflow and we set a free outflow boundary. In this case, we assume the two Riemann 
invariants are constant along the outgoing characteristics of the hyperbolic part (see 
[10]), therefore, we have a Neumann boundary condition for and 

5u"+1/2 


dH 

dn 


= 0 , 


dn 


= 0 . 


Preserving these conditions at the correction step, it gives a Neumann boundary 
condition for the pressure of type (47) 

dp 


dn 


= 0 . 


For an inflow given, the functional spaces will be defined by 


^ ^ ^out-fP ^ Qin' 


3. Numerical approximation. 

3.1. Discretization. This section is devoted to the numerical approximation and 
mainly for the correction step. Let us be given a subdivision of 17 with N vertices 
xi < X 2 < ■" < xn and we define the space step Aa:i_|_i /2 = Xi+i — Xi. We also note 
Ax^ = X^+i /2 - Xi_i /2 with Xi+i /2 = 

Prediction part. For the prediction step (20) i.e the hydrostatic part of the model, 
we use a finite volume scheme. We introduce the finite volume cells Ci centered at 
vertices Xi such that 17 = yJi=i^NCi. Then, the approximate solution A" at time 7" 

x{x,ndx, 

is solution of the numerical scheme 

\rn-\-l \rn ( -r-n rn \ i on 

J ^ 

where erf = and J- (resp. S) is a robust and efficient discretization of the 
conservative flux F{X) (resp. the source term 5'(X)). The time step is determined 
through a classical CFL condition. Many numerical fluxes and discretizations are 
available in the literature [8, 20, 29], we choose a kinetic based solver [5] coupled 
with the hydrostatic reconstruction technique [4] . 
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Correction part. Concerning the correction step (21)-(22), we consider the discrete 
problem corresponding to the mixed problem (40)-(41). We approach {Vo,Q) by the 
finite dimensional spaces {Voh,Qh) and we note 

N = dim(Vb/i), M = dim((5,i). 

We also denote by jy m basis functions of Voh and Qh re¬ 

spectively. The finite dimensional spaces will be specihed later on. We approximate 
(u,p) e {Vo,Q) by {uh,Ph) G {Voh,Qh) such that 

^hix) = J2[ Ph{x) ='^piMx)- 

z=l ' l=l 

Therefore, we consider the discrete problem: 


Find uji G Vb/i, Ph G Qh such that 


—a[\lh,Vh) + h{^rh,Ph) = 
bi^h,qh) = 

Let us introduce the mass matrix Mfj 


G Qhi 


At” ^ ^ 

0 , 


given by 


Vvft, G Vqh, 


(50) 

(51) 


Mh 



Hipiipjdx 


1 < i,i < N 


and the two matrices B*, B defined by 


B*" = i I Vswi4‘i)Pidx 


and we denote 


,B = - 




U = 


/ Ml \ 

UN 

Wi 


P = 


\wn/ 

Therefore, the problem (50)-(51) becomes 


B 


with 


0 


Ah = 


— At 


u 

p 


Mh 0 
0 Mh 


' diVsw{Po)(t>idx 

o 

/ Pi 
\PM/ 


l<l<Mp<j<N 


Assuming that Mh is invertible and eliminating the velocity U, we obtain the 
following equation 


BA7}B*P = -J—BU^+^/'^, 


(52) 


that is a discretization of the elliptic equation (44) of Sturm-Liouville type governing 
the pressure p. 
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We now take into consideration the boundary conditions in the more general 
problem (45)-(46). The velocity u is approximated by S 14, and the discrete 
problem is then written: 

Find {uh,Ph) e {Vh,Qh) such that 

-^a{uh,Vh) + b{vh,Ph) = -^a{ul'^^^^,Vh) + ^c{uh,Ph), Vv?, e 14, 

b{uh,qh) = 0 yqii&Qh- 

Considering the matrix = (c (pi, 4 ‘i))i<i<m i<i<N contains the bound¬ 

ary terms, the equation (52) becomes 

- C)P = 

This approach is suitable for the finite element approximation that is given in the 
next section. However, it implies to inverse a mass matrix Mfj that is not diagonal 
and depends on the water depth H. In practice, we use the mass lumping tech¬ 
nique introduced by Gresho ([22]) to avoid inverting the mass matrix in projection 
methods for Navier-Stokes incompressible system. 

3.2. Finite element Pi/Pg. In this part, the problem is solved by the mixed Hnite 
element approximation Pi/Pg (see [32]) on the domain fl = yjfL^Ki ( M = 1 with 

N the number of nodes), where the velocity is approximated by a continuous linear 
function and the pressure is approximated by a discontinuous piecewise constant 
function over each element 


w e 14 = {v„ e (C°(H))2 I g P?,V/ = 1,..., iv - 1}, 

and 

Ph&Qh = {qh I qhlxi ^ 'P’o , VI = 1 ,..., M - 1}. 

Using the discretization given in 3.1, we denote by iVi+ 1/2 the finite element [xi, cci+i], 
then the pressure is constant on the finite element IVi+i/ 2 . 

For the sake of clarity, in this situation, let ('('j-i- 1 / 2 ) be the basis functions 

for the pressure ph, and (</5i)i<j<jv the basis functions for the velocity u/i such that 



M 

Ph{x) = J2Pj + l/2<l^j+l/2{x). 

i=i 


We note (^ = H + 2zi, and assume () is approximated by a piecewise linear function 
Ch, namely Chix) = Ci^iix)■ We also note ^li+ 1/2 = = Xz+i /2 the 


constant gradient of C/t on the element 

If we denote p = {p, p)*, then the shallow water gradient operator is written 


XIsw{Ph) ■ Pi dx = 


_ f Hi{Pi+l /2 — Pi- 1 / 2 ) + Xi- 1/2 + ^‘^ 2 ^^^Xi+l /2 


— {IXXi+i/2Pi+l/2 + ^Xi-\l2Pi-\l2) 


Similarly, the shallow water divergence operator writes 
/ divs^ (u;i) 4 j+i /2 dx = Uj+i - Uj + Axj+i/ 2 {wj + Wj+i). 

J 

In one dimension, this approach corresponds to a staggered-grid finite-difference 
method where the velocity is computed at the nodes and the pressure is computed 
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at the middle nodes. The discretization we obtain corresponds exactly to the fi¬ 
nite difference scheme given in [1], and then, the properties established in [1] are 
conserved. 

3.3. Finite element Pi-iso-P 2 /Pi. For the one dimensional Pi-iso-P 2 /Pi, we con¬ 
sider two meshes ICh (tbe same as before) and lC 2 h with = [xi,Xi+i] and 

K 2 h,j = [x 2 j-i,X 2 j+i] the finite elements defined on the respective meshes K-h and 
IC 2 h such that ICh = C)^^Khh+i /2 ^ud IC 2 h = with N the total number 

of vertices of ICh and M = {N — l)/2 (assuming N odd), the number of vertices of 
lC 2 h- Therefore, the approximation spaces 14, and Qh are defined by 


Vh = 

{v,ec°(+ 2 | V,|^^_,GP 2 ,W = 1 , 

,A-l} 

Qh = 

{g, ePi,Vj = l,. 



Then, the velocity and the pressure are written 

M N 

Ph(x) = Uh{x) = 

j=l i=l 



Pi 



X2j — \ Xi~2j 


Figure 2. Representation of the basis functions. 


In figure 2, the dashed lines are the usual elementary basis functions of Pi on 
the mesh Kh, while the continuous lines are the basis functions on the mesh K 2 h- 
We can define the divergence operator, for all j = 1, M 


divsu,(u?,)0j dx = 


E 


divsu,(u/i)(/)j dx. 


IK 


We use a linear interpolation for Hipj^ and consider that Ati = Ax Vi = 1,..., iV 
for the sake of simplicity. We still approximate C by C,h defined before. 

The discrete shallow water divergence operator is computed for all nodes Xj of 
the mesh IC 2 h and therefore, denoting i = 2j — 1, it can be written, Vj = 1,M 


ti+2Ui+2 


[ divsu,(uii)(/),■ dx = ( -j 

( *-1/2 , i-l-l/2\ ( J-3/2\ / 

- j + Xi+i/2mij ' JUi- (^Xi-3/2m,_2j J u*_2 - (^X*+3/21 l 

/ i-1/2 , i-3/2\ ( i+1/2 , i+3/2 

\Xi-\l2m^-ij +Xi-3/2Ui,_iV [x^+l/2m,^l[J + Xz+3/2m,^i'j 

/ 1 — 1/9 i — 'iJ-‘CK/o\ 


■ m 


I / 2—3/2\ I 2+3/2 

+(’^i-2j + ni,+2j 'a'*+2i 


i-3/2^ 


Wi-l + 




i+3/2\ 


+ j 


with = + 

^ h, i + 1/2 


ipi(j)jdx. 


Ui+2 
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Similarly, the gradient shallow water operator is obtained for all the nodes Xi of 
the mesh K.h- However, we distinguish the gradient at the nodes of the elements 
K 2 h from the ones at the interior. In other words, for all the nodes Xi of the mesh 
lC 2 h, the gradient operator is defined by 


^swPh- 




H., 


= {Pj+i - Pj-i) 


4 

+Xi-l/2 




z , t — 1 / z 

Pj 


, / i+l/2 , i+1/2 \ 

+Xj+i/2 [rriij Pj + Pj+i j , 

n 2+1/2 2—1/2 2+1/2 

-2to, , Pj - TO, . Pj_i - TO, . Pj + I. 


^,3 


On the other hand, for all the nodes Xi such that i is even 


^sw Ph ■ 


'^sw Ph ■ P_ii^2j) 


H, 


-{Pj+i -P]) 

2 — 1/2 
X*-l/ 2 Wij 


i+1/2 


+ ( 

, ( i-1/2 , i+l/2\ 

+ (^X*-l/2TOiy+l + X*+l/2m,_^-+ JPl + l, 




- 1/2 




i+1/2 


\ / 1-1/2 , J+l/2\ 

) Pj - + ™*,1 + 1 ) Pj + i- 


With the discretization of the shallow water operators given below, we are able 
to validate the scheme for the two first order methods. Nevertheless, notice that in 
the following section, only the first order has been implemented. 


4. Validation with an analytical solution. In [1],[I2] some analytical solutions 
of the model (2)-(5) have been presented and they allow to validate the numerical 
method. We consider the propagation of a solitary wave without topography. This 
solution has the form 


H 

— Uq + a 

u 

= Co ^I - 

w 

acod 

= iff'’ 


X — Cot 

I 


seen 


P = 


u ,11 (x-Cot 

Hsech -;- sech 


aq 


with d, u.Hq gR, Hq > 0, a > 0 and cq 


^ / = 


dV 



The solitary wave is a particular case where dispersive contributions are coun¬ 
terbalanced by non linear effects so that the shape of the wave remains unchanged 
during the propagation. The propagation of the solitary wave has been simulated 
for the parameters a = OArrijHo = Ito, and d = Im over a domain of 45 to with 
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9000 nodes. At time t = 0, the solitary wave is positioned inside the domain. The 
results presented in figure 3 show the different fields, namely the elevation, the com¬ 
ponents of velocity and the total pressure at different times, and the comparison 
with the analytical solution at the last time. 




Figure 3. Propagation of the solitary wave at times 
1.00008 s, 1.9009 s, 3.9017 s and 5.9025 s. Comparison with 
analytical solution at time t = 5.9025 s. 

In the projection step, the greatest difficulty is to compute the pressure corre¬ 
sponding to the boundary conditions of the hyperbolic part (as seen in 2.3). The 
solution near the boundary has been confronted to the analytical solution. In the 
following result, we set a Neumann boundary condition on the non hydrostatic pres¬ 
sure with the parameters given below. As shown in the figure 4, the pressure is well 
estimated at the outflow boundary and allows the wave to leave the domain with 
a good behavior. The inflow boundary condition has been tested with this same 
test case and gives similar results. We are able to let the solitary wave enter in the 
domain with a good approximation of the elevation. 

The numerical simulations for the first order method are compared with the 
analytical solution and the error has been evaluated over different meshes of 
sizes from 603 nodes to 6495 nodes (see figure 5). With the parameters given 
above, it gives a convergence rate close to 1 for the two computations, i.e Pi/Pq 
and Pi-isoP2/Pi. 

Notice that the parameters set to validate the method lead to have a significant 
non hydrostatic pressure (see the figure 4) and then, the results show the ability of 
the method to preserve the solitary wave over the time. The numerical results have 
also been obtained for the Thacker’s test presented in [1], with the same rate as the 
Pi/Po method. 
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Figure 4. Non hydrostatic pressure profile at right boundary {x = 
45m) (a): t = 9.4044s (b): t = 9.8046s (c): t = 10.1048s. 



Figure 5. Convergence rate of the solitary wave solution for Pi/Pq 
and Pi-iso-P2/Pi scheme. 


5. Numerical results. 

5.1. Dam break problem. We next study the dispersive effect on the classical 
dam break problem, which is usually modeled by a Riemann problem providing 
a left state {H[^,ul) and a right state (HnjUn) on each side of the discontinuity 
Xd ([20]). However, our numerical dispersive model does not allow discontinuous 
solutions due to the variational spaces required for H (see also [12]), thus we provide 
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an initial data numerically close to the analytical one 

H{x,0) = {Hu + a) — atanh ^— — 

a = Hr-Hl. 

To evaluate the non hydrostatic effect, the different fields have been compared 
with the shallow water solution with the initial data: = 1.8 to, Hr = Ito, 

Ur = Ur = e = 10“^ TO, Xd = 300 TO over a domain of length 600 to with 

30000 nodes. In figure 6, the evolution of the state is shown at time t = 10 s and 
t = 45 s. The oscillations are due to the dispersive effects but the mean velocity 
does not change. These results are in adequation with the analysis proposed by 
Gavrilyuk in [28] for the Green-Naghdi model with the same configuration. 



Figure 6. The dam break problem, elevation H and velocity u at 
times t = 10 s and t = 45 s. 


5.2. Wet-dry interfaces. The ability to treat the wet/dry interfaces is crucial in 
geophysical problems, since geophysicists are interested in studying the behavior of 
the water-depth near the shorelines. This implies a water depth tending to zero 
at such boundaries. To treat the problem, we use the method introduced in [1], 
considering a minimum elevation H^. 

Therefore, we confront the method with a coastal bottom at the right boundary 
over a domain of 35 to with 3000 nodes. A wave is generated at the left boundary 
with an amplitude of 0.2 to and an initial water depth Hq = 1 to. In figure 7, the 
arrival of the wave at the coast is shown for times t = 7.91 s, 9.92 s and 10.42 s. 
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Figure 7. Propagation of a wave at a wet/dry interface. 

5.3. Comparison with experimental results. In this part, we confront the 
model with Dingemans experiments (detailed in [18, 17]) that consist in gener¬ 
ating a small amplitude wave at the left boundary of a channel with topography as 
described in figure 8. 



Figure 8 . Configuration of Dingemans’s test. 

At the left boundary, a wave is generated with a period T = 2.02 s and an 
amplitude of 0.02 m. A free outflow condition is set at the right boundary. The 
initial free surface is set to be rjQ = 0.4 m, and the measurement readings are saved 
at the following positions 10.5 m, 12.5 m, 13.5 m, 14.5 m ,15, 7 m and 17.3 to, placed 
at sensors 1 to 6 (fig. 8 ). In such a situation, the non hydrostatic effects have a 
significant impact on the water depth that cannot be represented by a hydrostatic 
model. These effects result mainly from the slope of the bathymetry, 10% in this 
case. In the figure 9, the simulation has been run with the hydrostatic model and 
the elevation has been compared with measures at the sensor 5. As one can see, 
the non-hydrostatic pressure has to be taken in consideration to estimate the real 
water depth variation. 

The numerical simulation with the non-hydrostatic model has been run with 
15000 nodes on a domain of 49 to over 25 s and the comparisons are illustrated for 
each sensor (fig. 10). The goal of this last result is also to highlight the ability of 
the model to capture dispersive effects for a geophysical flow with a non negligible 
pressure. 

5.4. Remark on iterative method. We recall that this formulation should allow 
to extend the method on two dimensional unstructured grids. However, it requires 
to inverse a system at each time iteration, which will become too costly in two 
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Figure 9. Comparison with hydrostatic model on sensor 5. 






dimensions. To anticipate the two dimensional problem, this method has been 
tested using different iterative methods like conjugate gradient and Uzawa methods. 
In hgure 11, we show a comparison of the computing time for the implementation 
Pi-iso-P2/Pi and Uzawa method. In one dimension, it is not relevant to use one of 
these methods, while it will be necessary for the two dimension model. 

6. Conclusion. In this paper, a variational formulation has been established for 
the one dimensional dispersive model introduced in [12]. The main idea is to give 
a new framework in which it will be possible to extend the scheme to the two 
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Figure 10. Comparison between measured and computed eleva¬ 
tions on Dingemans test for the first six sensors. 



Figure 11. Comparison of the computing time (CPU) for the 
Pi/Po scheme, the Pi-iso-P 2 /Pi and Uzawa method with Pi-iso- 
P 2 /P 1 solver with a tolerance 10“^. 


dimensional model. To this aim, the finite-element method has been presented 
with two approximation spaces. First, the Pi/Pq approximation has been done and 
we recover, as expected, the hnite difference scheme, together with the good results 
proved in [1]. Then, the Pi-iso-P 2 /Piapproximation has been studied to prepare the 
two dimensional problem. We have validated the method using several numerical 
tests and studying the dispersive effect on geophysical situations. 
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